
import cv2
import math
import numpy as np
import matplotlib.pyplot as plt
import random
from sklearn.cluster import KMeans


'''helper function to display images'''
def display_images(images,titles,row,col):
    fig = plt.figure(figsize = (20,20))
    for i in range(len(images)):
        fig.add_subplot(row,col,i + 1)
        #RGB
        if np.ndim(images[i]) == 3:
            plt.imshow(images[i])
        else:
            plt.imshow(images[i],cmap = 'gray')
        plt.title(titles[i])
        plt.axis("off")
        
        
''' function load a input image,flash image,trimap h for an image and display the images '''
def load_images(index):
    f0 = cv2.imread("./../data_3/" + str(index) + "/noflash.png",cv2.IMREAD_UNCHANGED)
    f1 = cv2.imread("./../data_3/" + str(index) + "/flash.png",cv2.IMREAD_UNCHANGED)
    tm = cv2.imread("./../data_3/" + str(index) + "/trimap.png",cv2.IMREAD_UNCHANGED)
    #print(im.shape,fl.shape,tr.shape)
    f0 = cv2.cvtColor(f0,cv2.COLOR_BGR2RGB)
    f1 = cv2.cvtColor(f1,cv2.COLOR_BGR2RGB)
    f2 = f1 - f0
    #print(f0.shape,f1.shape,f2.shape,tm.shape)
    display_images([f0,f1,f2,tm],["no flash","flash","only flash","trimap"],2,2)
    return f0,f1,f2,tm

''' helper function to divide the trimap into surely foreground and surely background and unknown region'''
def get_maps(trimap):
    bg = np.copy(trimap.astype(np.uint8))
    fg = np.copy(trimap.astype(np.uint8))
    un = np.copy(trimap.astype(np.uint8))
    bg[trimap != 0] = 0
    bg[trimap == 0] = 255
    fg[trimap != 255] = 0
    un[trimap == 125] = 255
    un[trimap != 125] = 0
    return bg,fg,un

''' helper to zero pad the image '''
def pad_image(img,pad_x,pad_y):
    #grayscale
    if(np.ndim(img) == 2):
        (l,b) = img.shape
        new_image = np.zeros((l + 2 * pad_y, b + 2 * pad_x),dtype = np.uint8)
        new_image[pad_y : pad_y + l,pad_x : pad_x + b] = img
        return new_image
    #RGB
    else:
        (l,b,h) = img.shape
        new_image = np.zeros((l + 2 * pad_y, b + 2 * pad_x , h),dtype = np.uint8)
        new_image[pad_y : pad_y + l,pad_x : pad_x + b,:] = img
        return new_image

    ''' create a gaussian kernel of size k with sigma '''
def create_gaussian(k,sigma):
    return np.dot(cv2.getGaussianKernel(k,sigma),cv2.getGaussianKernel(k,sigma).T)


''' helper function to get neighbourhood '''
''' assuming the height and width to be odd for convenience'''
''' to do handle exceptions for now assuming its correct'''
def get_neighbourhood(img,centre_y,centre_x,height ,width ):
    new_centre_y = centre_y + height // 2
    new_centre_x = centre_x + width // 2
    ''' grayscale'''
    #print(np.ndim(img))
    if(np.ndim(img) == 2):
        pad_x = width // 2
        pad_y = height // 2
        padded_img = pad_image(img,pad_x,pad_y)
        (l,b) = img.shape
        left = new_centre_x - width // 2
        right = new_centre_x + width // 2
        up = new_centre_y - height // 2
        down = new_centre_y + height // 2
        if(left < 0 or up < 0 or right >= padded_img.shape[1] or down >= padded_img.shape[1] ):
            ''' for debugging 
            print("Neighbourhood error")
            print("Neighbourhood")
            print(left,right,up,down)
            print("Shapes")
            print(padded_img.shape,img.shape,height,width)
            print("Centre : {0},{1}".format(centre_y,centre_x))
            print("New Centre : {0},{1}".format(new_centre_y,new_centre_x))
            print("Done")
            '''
            print("Invalid indexes for computing neighbourhood")

        return padded_img[up : down + 1,left : right + 1]
   
    elif np.ndim(img) == 3:
        
        (l,b,h) = img.shape
        pad_x = width // 2
        pad_y = height // 2
        padded_img = pad_image(img,pad_x,pad_y)
        left = new_centre_x - width // 2
        right = new_centre_x + width // 2
        up = new_centre_y - height // 2
        down = new_centre_y + height // 2
        if(left < 0 or up < 0 or right >= padded_img.shape[1] or down >= padded_img.shape[1] ):
            ''' for debugging 
            print("Neighbourhood error")
            print("Neighbourhood")
            print(left,right,up,down)
            print("Shapes")
            print(padded_img.shape,img.shape,height,width)
            print("Centre : {0},{1}".format(centre_y,centre_x))
            print("New Centre : {0},{1}".format(new_centre_y,new_centre_x))
            print("Done")
            '''
            print("Invalid indexes for computing neighbourhood")
        return padded_img[up : down + 1,left : right + 1,:]


f0,f1,f2,tm = load_images(1)


print(f0.shape)

(573, 943, 3)


#clusters the neighbourhood around a location and tries to compute the mean and convariance for each

def getMeanandCovariance(pixels,kernel,num_cluster = 5):
    #print("MEAN ")
    #print(pixels.shape)
    #shape of the pixels is (N,3) where N is the valid foreground/background pixels
    means = []
    Covs = []
    clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
    labels = np.array(clusters.predict(pixels))
    for i in range(num_cluster):
        group = (labels == i) #truth array of the pixels withing the cluster
        this_pixels = np.asarray(pixels[group]) #pixel values of the cluster
        kernel_bin = np.reshape(kernel[group], (kernel[group].size,1)) #reshape the weights to the this_pixel
        kernel_sqrt = np.sqrt(kernel_bin) # useful for covariance computing

        #compute the weighted mean
        #compute the mean for each channel individually
        meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),
                            np.dot(this_pixels[:,2], kernel_bin)])/np.sum(kernel_bin)
        #print("mean shape : {0} this_pixels shape : {1} kernel_bin.shape : {2}".format(meanbin.shape,this_pixels.shape,kernel_bin.shape))
        meanbin = np.reshape(meanbin, (3,))
        means.append(meanbin)

        diff = (this_pixels - meanbin)
        diff = np.multiply(kernel_sqrt, diff)
        covbin = np.eye(3)
        with np.errstate(divide='raise'): # raise errors, but dont stop execution when low number of clusters, and other cases
            try:
                covbin = np.dot(diff.T, diff)/(np.sum(kernel_bin) + 0.00001) + np.eye(3)*(10**(-5)) #adding small value to prevent singular matrices
            except Exception as e:
                print("Error")
                print(covbin.shape)
            #print("Maybe")
            #print(covbin)
        Covs.append(covbin)  
    means = np.array(means).T
    Covs = np.array(cv2.merge(Covs)) #merge 3 channels to 1
    #means have shape (3,K) K is the number of cluster
    #Covs have the shape (3,3,K) where K is the number of cluster
    #print(means.shape,Covs.shape,num_cluster)
    return means, Covs


'''
Solves the matting Equations Iteratively
Choose the best pair of (FG,BG) which maximizes the likelyhood of a given F,B,alpha
returns the value of F,B, alpha
'''
def Solver(f0mean, f0Cov, bgmean, bgCov, f1mean, f1Cov, fmean, lmean, CCov, initial_alpha, minLikelihoodDelta = 10**(-3),maxIterations = 50):

    #means have shape (3,K) K is the number of cluster
    #Covs have the shape (3,3,K) where K is the number of cluster
    #CMean is a [a,b,c]
    #CCov is a scaler
    #initial_Alpha is a float
    #Some global required constants
    z,num_cluster = f0mean.shape
    solved_f0 = np.zeros(3)
    solved_bg = np.zeros(3)
    solved_f1 = np.zeros(3)
    solved_alpha = 0
    I = np.eye(3)
    Z = np.zeros((3,3))
    nIterations = 0
    maxLikelihood = - np.inf
    
    for i in range(num_cluster):
        try:
            invCovF0 = np.linalg.inv(f0Cov[:,:,i]) # compute inverse of the FG convariance matrix
        except Exception as e: # maynot be singular
            print("The Foreground matrix is not singular,Skipping for now")
            continue

        cur_f0mean = f0mean[:, i] #mean of ith FG Cluster 1*3
        
        for j in range(num_cluster):
            try:
                invCovF1 = np.linalg.inv(f1Cov[:,:,j]) # compute inverse of the FG convariance matrix
            except Exception as e: # maynot be singular
                print("The Foreground Flash matrix is not singular,Skipping for now")
                continue

            cur_f1mean = f1mean[:, j] #mean of ith FG Cluster 1*3

            for k in range(num_cluster):
                try:
                    invCovB = np.linalg.inv(bgCov[:,:,k]) # compute inverse of the BG convariance matrix
                except Exception as e: # Sometimes its singular, so we dump the matrix and skip for now!
                    print("The Background matrix is not singular,Skipping for now")
                    continue


                cur_bgmean = bgmean[:, k] #mean of ith BG Cluster

                #Initialization for the iterative solver
                alpha = initial_alpha
                nIterations = 0
                prevLikelihood = -np.inf
                while True:
                    nIterations += 1

                    #SOLVE THE EQUATION DESCRIBED IN THE PAPER

                    #define A
                    a11 = invCovF0 + I*(alpha/CCov)**2  #3X3
                    a12 = I*(alpha*(1-alpha)*(CCov)**2) # 3X3
                    a13 = Z
                    a22 = invCovB + I*((1 - alpha)/CCov)**2 #3X3
                    a33 = invCovF1 + I*(alpha/CCov)**2
                    r1 = np.hstack((a11, a12, a13)) #create the top row 3*9
                    r2 = np.hstack((a12, a22, a13)) # create the middle row 3*9
                    r3 = np.hstack((a13, a13, a33)) #create the bottom row
                    A = np.vstack((r1, r2, r3)) #create the final matrix 6*6

                    #define b
                    b11 = np.dot(invCovF0, cur_f0mean) + fmean*(alpha/CCov**2)  # 3*1
                    b12 = np.dot(invCovB, cur_bgmean) + fmean*(1-alpha)/(CCov**2) # 3*1
                    b13 = np.dot(invCovF1, cur_f1mean) + lmean*(alpha/CCov**2) #3*2
                    b = np.concatenate((b11, b12, b13)).T #(6,)

                    #solve Ax = b
                    try:
                    
                        #print(np.sum(A))
                        #print(np.sum(b))
                        '''
                        print("printing A")
                        for i1 in range(A.shape[0]):for j2 in range(b.shape[0]):
                            print("index {0} : {1}".format(j2,b[j2]))
                            for j1 in range(A.shape[1]):
                                print("index {0},{1} : {2}".format(i1,j1,A[i1,j1]))
                        print("printing B")
                        for j2 in range(b.shape[0]):
                            print("index {0} : {1}".format(j2,b[j2]))
                        '''
                        x = np.linalg.solve(A, b)
                    except Exception as e: # Sometimes there is an issue with solving if A is non invertible. This did not occur after singular matrices were skipped, but is still there for safety.
                        print("ERROR SOLVING")
                        print(e)
                        break


                    #assign f0,f1  and bg that are solved
                    f0Col = x[0:3]
                    bgCol = x[3:6]
                    f1Col = x[6:9]

                    # SOLVE FOR ALPHA using estimated F, B, C       
                    #print("f0Col : {0} bgCol : {1} fmean : {2} lmean : {3} f1col : {4}".
                    #      format(f0Col.shape,bgCol.shape,fmean.shape,lmean.shape,f1Col.shape))
                    num = np.dot((f0Col - bgCol).T, (fmean.T - bgCol)) + np.dot(f1Col.T,lmean)
                    den = np.dot((f0Col - bgCol).T, (f0Col - bgCol)) + np.dot(f1Col.T,f0Col)
                    #print(num/den)
                    alpha = num / den
                    #clip Alpha
                    if alpha < 0:
                          alpha = 0
                    if alpha > 1:
                          alpha = 1


                    # Compute Likelyhood

                    #L[0] = -np.sum((Cmean.T - alpha*fgCol - (1-alpha)*bgCol)**2)/(colCov**2) # L(C|F,B, alpha)
                    L1 = -(np.linalg.norm((fmean.T - alpha*f0Col - (1-alpha)*bgCol),2)**2)/(CCov**2) # L(C|F,B, alpha)
                    L2 = -(np.linalg.norm((lmean.T - alpha*f1Col),2)**2)/(CCov**2) # L(I'|F,, alpha)

                    LF0 = -(np.dot(np.dot((f0Col - cur_f0mean.T).T, invCovF0), (f0Col - cur_f0mean.T).T)/2) #L(F)
                    LF1 = -(np.dot(np.dot((f1Col - cur_f1mean.T).T, invCovF1), (f1Col - cur_f1mean.T).T)/2) #L(F)
                    LB = -(np.dot(np.dot((bgCol - cur_bgmean.T).T, invCovB), (bgCol - cur_bgmean.T).T)/2) #L(B)
                    L = np.array([L1,L2,LF0,LF1,LB])

                    likelihood = np.sum(L)

                    if likelihood > maxLikelihood: # If best likelihood so far, use that
                        solved_alpha = alpha
                        maxLikelihood = likelihood
                        solved_f0 = f0Col
                        solved_f1 = f1Col
                        solved_bg = bgCol

                    # Stop solving if the solver saturates or if it exceeds the maximum number of iterations per cluster pair
                    if abs(prevLikelihood - likelihood) < minLikelihoodDelta or nIterations > maxIterations: # Stop solving
                        break

                    prevLikelihood = likelihood

    return solved_f0, solved_f1 ,solved_bg, solved_alpha


#the main function,has responsible to process the whole Joint Bayesian Matting Algorith
# f0 is the RGB no flash image
# f1 is the RGB flash image
# f2 is the RGB only flash image given by f2 = f1 - f0
# trimap is the trimap generated by our code.
# threshold is the number of known pixels in the neighbourhood for to proceed with solving 

def solve(f0,f1,f2,tm,thresh = 10):
    #print(f0.shape,f1.shape,f2.shape,tm.shape)
    # get the regions for each type
    bg,fg,unk = get_maps(tm)
    y,x = np.where(unk != 0) # locations where classification is unknown
    unsolved_loc = list(set(zip(y,x)))
    unsolved_loc = sorted(unsolved_loc,key = lambda x : x[0]) #top to bottom
    # this is essentially the alpha matrix
    # 1 for foreground pixels
    # 0 for background pixels
    # nan for unknowns.
    mask = np.copy(unk).astype(np.float64)
    
   
    
    mask[unk != 0] = np.nan #replace unknown regions with nan
    mask[fg != 0] = 1 #replace foreground region with 1
    
    
    
    #convert fg to RGB

    fg1 = np.zeros((f0.shape))
    fg1[:,:,0] = fg
    fg1[:,:,1] = fg
    fg1[:,:,2] = fg
    #convert bg to RGB for bitwise operations
    bg1 = np.zeros((f0.shape))
    bg1[:,:,0] = bg
    bg1[:,:,1] = bg
    bg1[:,:,2] = bg
    
    #convert fg to RGB for bitwise
    fg2 = np.zeros(f1.shape)
    fg2[:,:,0] = fg
    fg2[:,:,1] = fg
    fg2[:,:,2] = fg

    #generate surely background and surely foreground images
    f = cv2.bitwise_and(fg1.astype(np.uint8),f0.astype(np.uint8))
    b = cv2.bitwise_and(bg1.astype(np.uint8),f0.astype(np.uint8))
    l = cv2.bitwise_and(fg2.astype(np.uint8),f2.astype(np.uint8))
    
    #display_images([f,b,im],["foreground","background","original"],3,1)
    
    iteration = 0
    solved = []
    tosolve = np.copy(unsolved_loc)
    pass_num = 0
    pass_thresh = 5 # after this threshold,window size need to be increased
    window_size = 15 # neighbourhood is (19,19)
    sigma = 10
    print("Unsolve : {0}".format(len(unsolved_loc)))
    while(len(tosolve) > 0):
        pass_num += 1
        tosolve = []
        #create a list of locations to be solved
        for loc in unsolved_loc :
            #check it is not already solved
            if loc not in solved:
                tosolve.append(loc)
        if pass_num > pass_thresh and pass_num % 3:
            window_size += 2 #initally odd,remains odd.Window size increases when too many unknowns
        print("PASS NUMBER : {0} , WINDOW SIZE : {1}".format(pass_num,window_size))
        for loc in tosolve:
            gaussian_kernel = create_gaussian(window_size,sigma) #create kernel
            alpha = get_neighbourhood(mask,loc[0],loc[1],window_size,window_size)
            
            #find  fg pixels for clustering non flash
            kernel_f0 = np.multiply(np.square(alpha),gaussian_kernel)
            f0_window = get_neighbourhood(f,loc[0],loc[1],window_size,window_size)
            known = np.nan_to_num(kernel_f0) > 0 #return a boolean array with true values for non-zero and non-nan values
            f0_pixel = f0_window[known] #extract the selected pixel location
            kernel_f0 = kernel_f0[known] # extract the selected pixel location
            
            #find  fg pixels for clustering flash
            kernel_f1 = np.multiply(np.square(alpha),gaussian_kernel)
            f1_window = get_neighbourhood(l,loc[0],loc[1],window_size,window_size)
            known = np.nan_to_num(kernel_f1) > 0 #return a boolean array with true values for non-zero and non-nan values
            f1_pixel = f1_window[known] #extract the selected pixel location
            kernel_f1 = kernel_f1[known] # extract the selected pixel location
            
            #find  bg pixels for clustering
            kernel_b = np.multiply(np.square(1- alpha),gaussian_kernel)
            bg_window = get_neighbourhood(b,loc[0],loc[1],window_size,window_size)
            known = np.nan_to_num(kernel_b) > 0 #return a boolean array with true values for non-zero and non-nan values
            bg_pixel = bg_window[known] #extract the selected pixel location
            kernel_b = kernel_b[known] # extract the selected pixel location
            
          
            
            #check sufficient data to solve the proble
            if(len(f0_pixel) >= thresh and len(bg_pixel) >= thresh) and len(f1_pixel):
                    

                f0mean,f0Cov = getMeanandCovariance(f0_pixel,kernel_f0,num_cluster = 3)
                f1mean,f1Cov = getMeanandCovariance(f1_pixel,kernel_f1,num_cluster = 3)
                bgmean,bgCov = getMeanandCovariance(bg_pixel,kernel_b,num_cluster = 3)
                
                #observed_c
                fmean = f0[loc[0],loc[1],:]
                lmean = f1[loc[0],loc[1],:]
                CCov = 3 #tunable parameter
                #print(Cmean,CCov)
                
                #take intial guess as mean
                initial_alpha = np.nanmean(alpha)
                
                #print(initial_alpha)
                
                
                #iteratively solve!
                final_f0, final_f1, final_bg , final_alpha = Solver(f0mean, f0Cov, bgmean, bgCov, f1mean, f1Cov, fmean, lmean, CCov, initial_alpha, maxIterations = 50)


                #update the estimated values
                
                #print(final_fg,final_bg,final_alpha)
                mask[loc[0],loc[1]] = final_alpha
                f[loc[0],loc[1],:] = np.array(final_f0)
                b[loc[0],loc[1],:] = np.array(final_bg)
                l[loc[0],loc[1],:] = np.array(final_f1)
                solved.append(loc)


                #print occasional status updates
                if iteration % 50 == 0:
                    print("PASS : ", pass_num)
                    print("WINDOW SIZE = ", window_size)
                    print('SOLVED LOCATIONS = ', len(solved))
                    print('REMAINING LOCATIONS = ', len(tosolve) - len(solved))
                    #sad = np.sum(np.absolute((cv2.cvtColor(gt, cv2.COLOR_BGR2GRAY)/255).astype(np.float) - np.nan_to_num(mask)))
                    #print('SAD = ', sad)
                    print('---------------------------------------------')
                    final_alpha = cv2.cvtColor((mask*255).astype(np.uint8), cv2.COLOR_GRAY2BGR) #3 c
                    #showMultiImages((fg, bg, gt, alpha3channel), 'CURRENT')


                iteration += 1
    #sad = np.sum(np.absolute((cv2.cvtColor(gt, cv2.COLOR_BGR2GRAY )/255).astype(np.float) - np.nan_to_num(mask)))
    print('SOLVER COMPLETED IN ', iteration, 'ITERATIONS')
    print('---------------------------------------------')
    print('---------------------------------------------')
    print('---------------------------------------------')


    return f, b, l,mask


def scale(im,sc):
    return cv2.resize(np.copy(im),(0,0),fx = sc,fy = sc)


def run(index):
    import time
    t1 = time.time()
    f0,f1,f2,tm = load_images(index)
    f0 = scale(np.copy(f0),0.5)
    f1 = scale(np.copy(f1),0.5)
    f2 = scale(np.copy(f2),0.5)
    tm = scale(np.copy(tm),0.5)
    f,b,l,mask = solve(f0,f1,f2,tm,thresh = 10)
    t2 = time.time()
    print("Time taken : {0}".format(t2 -t1))
    return f,b,l,mask


f,b,l,mask = run(1)
display_images([f,b,l],["non flash foreground","background","flash foregroun"],3,1)

Unsolve : 6065
PASS NUMBER : 1 , WINDOW SIZE : 15
PASS :  1
WINDOW SIZE =  15
SOLVED LOCATIONS =  1
REMAINING LOCATIONS =  6064
---------------------------------------------
PASS :  1
WINDOW SIZE =  15
SOLVED LOCATIONS =  51
REMAINING LOCATIONS =  6014
---------------------------------------------

<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (4) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),

PASS :  1
WINDOW SIZE =  15
SOLVED LOCATIONS =  101
REMAINING LOCATIONS =  5964
---------------------------------------------

<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (4) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),
<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (4) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),
<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (1) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),
<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (1) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),
<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (2) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),
<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (1) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),
<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (1) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),
<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (1) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),
<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (1) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),
<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (1) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),
<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (1) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),

PASS :  1
WINDOW SIZE =  15
SOLVED LOCATIONS =  151
REMAINING LOCATIONS =  5914
---------------------------------------------

<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (1) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),
<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (4) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),
<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (3) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),

PASS :  1
WINDOW SIZE =  15
SOLVED LOCATIONS =  201
REMAINING LOCATIONS =  5864
---------------------------------------------

<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (2) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),

PASS :  1
WINDOW SIZE =  15
SOLVED LOCATIONS =  251
REMAINING LOCATIONS =  5814
---------------------------------------------

<ipython-input-27-ede0654d1681>:9: ConvergenceWarning: Number of distinct clusters (2) found smaller than n_clusters (5). Possibly due to duplicate points in X.
  clusters = KMeans(n_init = 5, n_clusters = num_cluster , random_state = 0).fit(pixels)
<ipython-input-27-ede0654d1681>:19: RuntimeWarning: invalid value encountered in true_divide
  meanbin = np.array([np.dot(this_pixels[:,0], kernel_bin),np.dot(this_pixels[:,1], kernel_bin),


np.unique(mask)

array([0.        , 0.78136151, 0.78361439, 0.81567336, 0.81842579,
       0.82816523, 0.82930855, 0.8370937 , 0.83833196, 0.84654254,
       0.8516133 , 0.85555014, 0.86102867, 0.86269305, 0.86407806,
       0.89236169, 0.90015179, 0.90365709, 0.90366519, 0.91595168,
       0.93285823, 0.93549314, 0.93552378, 0.93764044, 0.95073973,
       0.9515216 , 0.9572191 , 0.97248966, 0.97313811, 0.97555714,
       0.99097726, 0.99153465, 1.        ])


f0 = scale(np.copy(f0),0.5)
plt.imshow(f0)

<matplotlib.image.AxesImage at 0x7f5f01646ee0>


f0 = scale(np.copy(f0),0.5)
final_image = np.copy(f0)
for i in range(mask.shape[0]):
    for j in range(mask.shape[1]):
        if(mask[i,j] <= 0.75):
            final_image[i,j,:] = 0
plt.imshow(final_image)

<matplotlib.image.AxesImage at 0x7f5f01e653a0>